(* Find solution of f(x). Plot and sort out the number of solutions with bisection *)
Remove["Global`*"]
f[x_]:=x^3+x-1;
Plot[f[x],{x,0.5,0.8}]
{a,b}={0,1};
epsilon = 0.001;
While[b-a>epsilon,
c=0.5*(a+b);
Print[{a,c,b}];
If[f[a]*f[c]<0,b=c,a=c]
]
NSolve[f[x]==0,x][[3]]